
# Save time-average (CAGR) growth rates for visualization
t.plot = rownames(global_growth_percap)
cagr.growth.rates = apply(global_growth_percap, 2, function(x) cumsum(x)/1:length(t.plot))

## Graph growth distribution - Figure 1
pdf(figure.dir%&%'/Average_Growth_Rate_Distribution_'%&%today()%&%'.pdf',
    family='serif', width=10, height=7) #width=8.5,height=11)
plot(apply(cagr.growth.rates, 1, mean) ~ t.plot, type='l', lwd=2, ylim=c(-0.02,0.08),
     main='', ylab='Growth Rate (Time Average from 2020 to Year)', xlab='Year', 
     cex.lab=1.5, cex.axis=1.5)
grid(col='gray70')

axis(1, cex.axis=1.5, at=t.plot[1]); abline(v=t.plot[1], lty=3, col='gray70')
alpha.temp  = 0.4
# set.seed(2)
# g.range = apply(cagr.growth.rates[,sample(1e4, 1e3, replace=T)], 1, quantile, probs=c(0.005, 0.995,0.25, 0.75, 0.1, 0.9, 0.025, 0.975, 0.05, 0.95, 0.01, 0.99))
g.range = apply(cagr.growth.rates, 1, quantile, probs=c(0.005, 0.995,0.25, 0.75, 0.1, 0.9, 0.025, 0.975, 0.05, 0.95, 0.01, 0.99, 0.5))
dimnames(g.range)
alpha.temp  = 0.4
polygon(c(colnames(g.range), rev(colnames(g.range))), c(g.range['25%',], rev(g.range['75%',])), col=rgb(0,1,0, alpha.temp), border=NA )
polygon(c(colnames(g.range), rev(colnames(g.range))), c(g.range['75%',], rev(g.range['90%',])), col=rgb(0.25,0.5,0, alpha.temp), border=NA )
polygon(c(colnames(g.range), rev(colnames(g.range))), c(g.range['10%',], rev(g.range['25%',])), col=rgb(0.25,0.5,0, alpha.temp), border=NA )
polygon(c(colnames(g.range), rev(colnames(g.range))), c(g.range['5%',], rev(g.range['10%',])), col=rgb(0.5,0.25,0, alpha.temp), border=NA )
polygon(c(colnames(g.range), rev(colnames(g.range))), c(g.range['90%',], rev(g.range['95%',])), col=rgb(0.5,0.25,0, alpha.temp), border=NA )
polygon(c(colnames(g.range), rev(colnames(g.range))), c(g.range['95%',], rev(g.range['97.5%',])), col=rgb(1,0,0, alpha.temp*2/3), border=NA )
polygon(c(colnames(g.range), rev(colnames(g.range))), c(g.range['2.5%',], rev(g.range['5%',])), col=rgb(1,0,0, alpha.temp*2/3), border=NA )
# polygon(c(colnames(g.range), rev(colnames(g.range))), c(g.range['1%',], rev(g.range['99%',])), col=rgb(1,0,0, alpha.temp*2/3), border=NA )

abline(h=0, v=2021)
lg = legend('topright', legend=c('95% Range','90% Range', '80% Range', '50% Range'), #lty=c(1,NA,NA,NA,NA),
            fill=c(rgb(1,0,0, 0.5*2/3), rgb(0.5,0.25,0, 0.5), rgb(0.25,0.5,0, 0.5), rgb(0,1,0, 0.5)), bg='white')
lg = legend(x=lg$rect$left + 7, y=lg$rect$top-lg$rect$h, legend=c('Mean'), lty=1:2, lwd=2, bg='white')
dev.off()

